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ABSTRACT 

We have constructed all-sky Compton parameters maps (y-maps) of the thermal Sunyaev-Zeldovich (tSZ) effect by applying specifically tailored 
component separation algorithms to the 30 to 857 GHz frequency channel maps from the Planck satellite survey. These reconstructed y-maps are 
delivered as part of the Planck 2015 release. The y-maps are characterized in terms of noise properties and residual foreground contamination, 
mainly thermal dust emission at large angular scales and CIB and extragalactic point sources at small angular scales. Specific masks are defined 
to minimize foreground residuals and systematics. Using these masks, we compute the y-map angular power spectrum and higher order statistics. 
From these we conclude that the y-map is dominated by tSZ signal in the multipole range, 20 < £ < 600. We compare the measured tSZ power 
spectrum and higher order statistics to various physically motivated models and discuss the implications of our results in terms of cluster physics 
and cosmology. 

Key words, cosmological parameters - large-scale structure of Universe - Galaxies: clusters: general 
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1. Introduction 

This paper, one of a set associated with the 2015 release of data 
from the Planck^ mission, describes the Planck Compton param¬ 
eter (y) map, which is part of the Planck 2015 data release. 

The thermal Sunyaev-Zeldovich (tSZ) effect (Sunyaev & 
Zeldovich, 1972) is produced by the inverse Compton scatter¬ 
ing of cosmic microwave background (CMB) photons by hot 
electrons along the line of sight and in particular in clusters of 
galaxies. The tSZ effect has proved to be a major tool to study 
the physics of clusters of galaxies as well as structure forma¬ 
tion in the Universe. Catalogues of clusters of galaxies selected 
via the tSZ effect have become available in the last few years, 
including for example those from the Planck satellite (Planck 
Collaboration VIII, 2011; Planck Collaboration XXIX, 2014), 
the Atacama Cosmology Telescope (ACT, Hasselfield et al., 

2013) and the South Pole Telescope (SPT, Reichardt et al., 2013; 
Bleem et al., 2014). These catalogues and their associated sky 
surveys have been used to study the physics of clusters of galax¬ 
ies (Planck Collaboration XII, 2011; Planck Collaboration XI, 
2011; Planck Collaboration X, 201 1) and their cosmological im¬ 
plications (Planck Collaboration XX, 2014; Benson et al., 2013; 
Das et al., 2013; Wilson et al., 2012; Mak & Pierpaoli, 2012). 

The study of cluster number counts and their evolution with 
redshift using tSZ selected catalogues is an important cosmo¬ 
logical test (Carlstrom et al., 2002; Dunkley et al., 2013; Benson 
et al., 2013; Planck Collaboration XX, 2014; Hou et al., 2014). 
This study is well complemented by the measurement of the tSZ 
effect power spectrum (Komatsu & Seljak, 2002), for which no 
explicit measurement of cluster masses is required. Low mass 
clusters, thus fainter in tSZ, which can not be detected individu¬ 
ally, also contribute statistically to the measured signal (Battaglia 
et al., 2010; Shaw et al., 2010). Another complementary ap¬ 
proach (as pointed out in Rubino-Martm & Sunyaev, 2003) is to 
study the higher order statistics of the tSZ signal, and in partic¬ 
ular the skewness or, equivalently, the bispectrum. The bispec¬ 
trum of the tSZ effect signal is dominated by massive clusters 
at intermediate redshifts (Bhattacharya et al., 2012), for which 
high-precision X-ray observations exist. This contrasts with the 
power spectrum, for which most of the signal comes from the 
lower mass and higher redshift groups and clusters (e.g., Trac 
et al., 2011). Thus, theoretical uncertainties in the tSZ bispec¬ 
trum are expected to be significantly smaller than those in the es¬ 
timation of the tSZ power spectrum. Therefore, combined mea¬ 
surements of the power spectrum and the bispectrum can be used 
to distinguish the contribution to the power spectrum from differ¬ 
ent cluster masses and redshift ranges. However, contamination 
from point sources (Rubino-Martm & Sunyaev, 2003; Taburet 
et al., 2010) and other foregrounds (Planck Collaboration XX, 

2014) needs to be dealt with carefully. 

As shown in the Planck 2013 results (Planck Collaboration 
XX, 2014), the all-sky coverage and unprecedented wide 
frequency range of Planck allowed us to produce all-sky tSZ 
Compton parameter maps, also referred to as y-maps. From 
this map an accurate measurement of the tSZ power spectrum 
at intermediate and large angular scales, for which the tSZ 
fiuctuations are almost insensitive to the cluster core physics. 


^ Planck (http://www.esa.int/Planck) is a project of the 
European Space Agency (ESA) with instruments provided by two sci¬ 
entific consortia funded by ESA member states and led by Principal 
Investigators from Prance and Italy, telescope reflectors provided 
through a collaboration between ESA and a scientific consortium led 
and funded by Denmark, and additional contributions from NASA 
(USA). 



Fig. 1: Window functions corresponding to the spectral localisa¬ 
tion of the NILC and MILCA algorithms. For NILC (black) there 
are 10 Gaussian window functions defining 10 needlet scales. 
MILCA (red) uses 11 Gaussian overlapping windows. 


can be obtained. Furthermore, the expected non-Gaussianity 
properties of the tSZ signal can be studied using higher order 
statistical estimators, such as the skewness and the bispectrum. 
The y-map can be also used to extract the tSZ signal on regions 
centered at cluster positions and in particular to perform stack¬ 
ing analysis. Finally, this map may also be cross correlated with 
other cosmological probes (see for example Ma et al., 2014; Van 
Waerbeke et al., 2014; Hill & Spergel, 2014) as a consistency 
test. 

In this paper we construct revised tSZ all-sky maps from the 
individual Planck frequency maps. With respect to the Planck 
Compton parameter map in the Planck 2013 results (Planck 
Collaboration XX, 2014) we have extended the analysis to the 
full mission data set and performed an in depth characterization 
of its statistical properties. This extended analysis allows us to 
deliver this map as part of the Planck 2015 data release. The pa¬ 
per is structured as follows. Section 2 describes the Planck data 
used to compute the tSZ all-sky maps and the simulations used 
to characterize it. Section 3 presents the reconstruction of the 
Planck all-sky Compton parameter map. In Sect. 4 we discuss 
the validation of the reconstructed y-map in the pixel domain 
including signal and noise characterization. Section 5 describes 
the power spectrum analysis. Cross-checks using higher order 
statistics are presented in Sect. 6. The cosmological interpreta¬ 
tion of the results is discussed in Sect. 7, and we finally present 
our conclusions in Sect. 8. 

2. Data and simulations 

2.1. The Planck data 

This paper is based on the Planck's full mission, corresponding 
to five and eight full-sky surveys for the High Frequency 
Instrument (HFI) and Low Frequency Instrument (LFI) data 
respectively. 

The Planck channel maps are provided in HEALPix (Gorski 
et al., 2005) pixelization scheme at Vside = 2048. A noise map 
is associated with each channel map. This map is obtained from 
the half difference of maps made from the first and second half 
of each stable pointing period (also called ring). In the half- 
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NILCtSZmap 



MILCA tSZ map 



Fig. 2: Reconstructed Planck all-sky Compton parameter maps for NILC (top) and MILCA (bottom) in orthographic projections. 
The apparent difference in contrast observed between the NILC and MILCA maps comes from differences in the residual foreground 
contamination and from the differences in the filtering applied for display purposes to the original Compton parameter maps. For the 
MILCA method filtering out low multipoles reduces significantly the level of foreground emission in the final y-map. The wavelet 
basis used in the NILC method was tailored for tSZ extraction. For details see Planck Collaboration XXII (2015). 


difference maps the astrophysical emission cancels out, which 
makes them a good representation of the statistical instrumental 
noise. These half-difference maps are used to estimate the 
noise in the final Compton parameter map. In addition, survey 
maps, which are also available for each channel, will be used to 
estimate possible residual systematic effects in the y-map. 


For the purpose of this paper we approximate the Planck ef¬ 
fective beams by circular Gaussians, the FWHM estimates of 
which are given in Table 1 for each frequency channel. 

2.2. Simulations 

We also use simulated Planck frequency maps obtained from 
the Full Focal Plane (FFP) simulations (Planck Collaboration 
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NILCtSZmap 


MLCAtSZmap 



Fig. 3: A small region of the reconstructed Planck all-sky Compton parameter maps for NILC (left) and MILCA (right) at interme¬ 
diate Galactic latitudes in the southern sky centered at (0, -45) degrees in Galactic coordinates. 


Table 1: Conversion factors for tSZ Compton parameter y to 
CMB temperature units and the FWHM of the beam of the 
Planck channel maps. 


Frequency Tcmb g(v) FWHM 
[GHz] [Kcmb] [arcmin] 


100 . -4.031 9.66 

143 . -2.785 7.27 

217 . 0.187 5.01 

353 . 6.205 4.86 

545 . 14.455 4.84 

857 . 26.335 4.63 


XII, 2015), which are described in the Planck Explanatory 
Supplement (Planck Collaboration ES, 2013). 

These simulated maps include the most relevant sky compo¬ 
nents at microwave and millimetre frequencies, based on fore¬ 
grounds from the Planck Sky Model (PSM, Delabrouille et al., 
2013): CMB, thermal SZ effect, diffuse Galactic emissions (syn¬ 
chrotron, free-free, thermal and spinning dust and CO), radio 
and infrared point sources, and the clustered CIB. The tSZ sig¬ 
nal was constructed using hydrodynamical simulations of clus¬ 
ters of galaxies up to redshift 0.3. For higher redshifts pressure 
profile-based simulations of individual clusters of galaxies ran¬ 
domly drawn on the sky have been added. The noise in the maps 
was obtained from realizations of Gaussian random noise in the 
time domain and therefore accounts for noise inhomogeneities 
in the maps. 


The simulation set also includes Monte Carlo noise-only re¬ 
alizations for each Planck channel map. These will also be used 
to estimate the noise properties in the final y-map. 

3. Reconstruction of the all-sky tSZ maps 

The thermal SZ Compton parameter (Sunyaev & Zeldovich, 
1972) in a given direction, fi, can be expressed as 

y{n)= f ne^^crxdi', (1) 

where is the Boltzmann contsant, me the electron mass, (Tt the 
Thomson cross-section, d^- the distance along the line of sight, n, 
and and Tq are the electron number density and temperature. 

In CMB temperature units the tSZ effect contribution to the 
Planck maps for a given observational frequency v is given by 

AT 

Tf; -=^(v)y. (2) 

7 CMB 

Neglecting relativistic corrections, g(y) = v coth(v/2) -4, with 
X = hvKkBTcMB)- Table 1 presents the conversion factors for 
Compton parameter to CMB temperature, Kcmb, for each fre¬ 
quency channel, after integrating over the instrumental bandpass 
(Planck Collaboration IX, 2014). 

3.1. Reconstruction methods 

The tSZ effect signal in the Planck frequency maps is subdom¬ 
inant with respect to the CMB and other foreground emissions. 
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Fig. 4: In and cross scan contributions in the NILC (top line) and MILCA (bottom line) y-maps. From left to right we present the 
original y-maps, and their in and cross scan contributions for a small region at intermediate Galactic latitudes in the southern sky 
centered at (0, -45) degrees in Galactic coordinates. 


By contrast to the CMB emission, the tSZ effect from galaxy 
clusters is spatially localized and leads to a highly non-Gaussian 
signal. CMB-oriented component-separation methods (Planck 
Collaboration XII, 2014) are not adequate to recover the tSZ 
signal. We therefore use specifically tailored component sepa¬ 
ration algorithms to reconstruct the tSZ signal from the Planck 
frequency channel maps as in Planck Collaboration XX (2014). 
These algorithms rely on the spatial localization of the differ¬ 
ent astrophysical components and on their spectral diversity to 
separate them. As in Planck Collaboration XX (2014) we con¬ 
sider here the MILCA (Modified Internal Linear Combination 
Algorithm, Hurier et al., 2013) and NILC (Needlet Independent 
Linear Combination, Remazeilles et al., 2011) methods. Both 
are based on the Internal Linear Combination (ILC) approach 
that searches for the linear combination of the input maps that 
minimizes the variance of the final reconstructed map under the 
constraint of offering unit gain to the component of interest (here 
the tSZ effect, whose frequency dependence is known). Both al¬ 
gorithms have been extensively tested on simulated Planck data. 

For both methods, the Planck HFI maps from 100 to 
857 GHz, convolved to a common resolution of 10', are used. 
In the case of NILC we also use the LFI data at large angular 
scales {i < 300). Similarly, for both methods the 857 GHz map, 
which traces the thermal dust emission on large angular scales, is 


only used for multipoles i < 300 to minimise residuals from IR 
point sources and clustered Cosmic Infrared Bacground (CIB) 
emission in the final y-maps. 

3 . 1 . 1 . NILC 

In the multi-component extensions of NILC (Delabrouille et al., 
2009; Remazeilles et al., 201 1), initially developed to extract the 
CMB, the weights for component separation (i.e., covariances) 
are computed independently in domains of a needlet decomposi¬ 
tion (in the spherical wavelet frame). The needlet decomposition 
provides localization of the ILC filters both in pixel and in mul¬ 
tipole space, allowing us to deal with local contamination con¬ 
ditions varying both in position and in scale. We imposed con¬ 
straints to remove the CMB contamination and preserve the tSZ 
effect. To avoid strong foreground effects, part of the Galactic 
plane (corresponding to about 2% of the sky) was masked be¬ 
fore applying NILC to the Planck frequency maps. 

The localisation in multipole space is achieved by using ten 
Gaussian window functions {/t;(^)}i<;<io as bandpass filters^, so- 
called needlet bands, allowing for smooth localisation in i (see 

^ Note that in Planck Collaboration XXI (2014) cosine window func¬ 
tions were used. 
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Fig. 1). NILC performs a weighted linear combination of the 
bandpass filtered Planck maps for each needlet scale indepen¬ 
dently. The localisation in the spatial domain is achieved by 
defining scale-dependent zones over the sky on which the covari¬ 
ance matrices and ILC weights are computed. More precisely, 
the pixel domain on which the covariance is computed, is de¬ 
fined from the smoothing of the product of the relevant needlet 
maps with a symmetric Gaussian window in pixel space whose 
FWHM depends on the needlet scale considered. This avoids 
artificial discontinuities at pixel edges. The localisation reduces 
the number of modes on which the statistics is computed, this 
may be responsible for an ‘TLC bias” due to chance correlations 
between SZ and foregrounds (Delabrouille et al., 2009). At the 
coarsest scale in particular (first needlet band), the area of the 
spatial localisation must be large enough to counterbalance the 
lack of modes in multipole space due to spectral localisation. In 
practice, the zones for spatial localisation are not pre-defined but 
their area is automatically adjusted to the needlet scale consid¬ 
ered . The ILC bias Z^ilc is related to the number of channels A^ch 
and to the number of modes Delabrouille et al. (2009) as 


^ILC = -O' 


2 

sz 


A^ch-1 

N 


( 3 ) 


masked in both maps, leaving 67% of the sky. Residual Galactic 
contamination is also visible as diffuse positive structures in the 
MILCA y-map. We can also observe a granular structure in the 
NILC y-map that corresponds to an excess of noise at large an¬ 
gular scales as discussed in Sect. 4.2. 

Weaker and more compact clusters are visible in the zoomed 
region of the Southern cap, shown in Fig. 3. Strong Galactic and 
extragalactic radio sources show up as negative bright spots on 
the maps and were masked prior to any scientific analysis, as 
discussed below in Sect. 4.4.1. We can again observe residual 
Galactic contamination around the edges of the masked area, 
which is more important for MILCA. Finally, we note in the 
NILC and MILCA y-maps the presence of systematic residu¬ 
als along the scanning direction that show up as stripes. These 
are the consequence of stripes in the original Planck maps (see 
Planck Collaboration VIII, 2015, for a detailed discussion). We 
discuss the level of residual stripes in Sect. 4.1. 

In addition to the full Compton parameter maps, we also pro¬ 
duce the so-called first and last (F and L hereafter) Compton pa¬ 
rameter maps from the first and second halves of the survey rings 
(i.e., pointing periods). These maps are used for power spectrum 
estimation in Sect. 5 as well as to estimate the noise properties 
in the y-maps (see Sect. 4.2). 


This offers the possibility to adapt Nm and Vch in order to con¬ 
trol the ILC bias to a set threshold, as discussed in Remazeilles 
et al. (2013): given both the number of channels and the num¬ 
ber of spectral modes in ^-space at the needlet scale considered, 
our NILC algorithm consistently computes the number of spatial 
modes (similarly, the zone area for spatial localisation) required 
to control the ILC bias. 

3 . 1 . 2 . MILCA 

MILCA (Hurier et al., 2013), when applied to the extraction of 
the tSZ signal, also uses two spectral constraints: preservation 
of the tSZ signal (the tSZ spectral signature discussed above is 
assumed); and removal of the CMB contamination in the final 
SZ map, making use of the well known spectrum of the CMB. 
In addition, to compute the weights of the linear combination, 
we have used the extra degrees of freedom in the linear system 
to minimize residuals from other components (two degrees of 
freedom) and from the noise (two additional degrees). The noise 
covariance matrix was estimated from the half-difference maps 
described in Section 2.1 . To improve the efficiency of the MILCA 
algorithm, weights are allowed to vary as a function of multipole 
i, and are computed independently on different sky regions. We 
have used 11 filters in £ space as shown in Figure 1 . These filters 
have an overall transmission of one, except for ^ < 8. For these 
large angular scales we have used a Gaussian filter to reduce 
foreground contamination. To ensure sufficient spatial localiza¬ 
tion for each required resolution the size of the independent sky 
regions was adapted to the multipole range. We used a minimum 
of 12 regions at low resolution and a maximum of 3072 regions 
at high resolution. 

3.2. Reconstructed Compton parameter map 

Figure 2 shows the reconstructed Planck all-sky Compton pa¬ 
rameter map for NILC (top panel) and MILCA (bottom panel). 
For display purposes, the maps are filtered using the procedure 
described in the first paragraph of Sect. 6.1. Clusters appear as 
positive sources: the Coma cluster and Virgo supercluster are 
clearly visible near the north Galactic pole. The Galactic plane is 


3.3. Comparison to other Compton parameter maps in 
the literature 

Van Waerbeke et al. (2014) and Hill & Spergel (2014) have used 
the Planck nominal data to reconstruct a Compton parameter 
map over a large fraction of the sky. As for the present work, 
they use an ILC approach imposing spectral constraints to have 
unitary gain for the SZ component and null contribution from 
CMB, but they do not use the spectral and spatial separation 
discussed above. Van Waerbeke et al. (2014) consider only the 
four HFI channels between 100 and 353 GHz and force a null 
contribution from dust emission. On their side Hill & Spergel 
(2014) also include the 545 GHz map, while using the 857 GHz 
channel only as a template for dust emission: a fiux cut is im¬ 
posed in order to build a mask that keeps only the 30% of the 
sky. They also apply a point source mask (radio and IR) be¬ 
fore computing the ILC coefficients reducing the final sky frac¬ 
tion to about the 25%. Furthermore, Van Waerbeke et al. (2014) 
and Hill & Spergel (2014) aim mainly at studying the cross¬ 
correlation with the gravitational lensing mass map from the 
Canada France Hawaii Telescope Lensing Survey (CFHTLenS) 
and the publicly-released Planck CMB lensing potential map, 
respectively. Their y-map reconstruction methods are tailored to 
fulfill these objectives and as a consequence they present larger 
overall foreground contamination, in the case of Hill & Spergel 
(2014) also a significantly smaller sky fraction. 

4. Pixel space analysis 

4 . 1. Stripes in the y-map 

As discussed in Sect. 3.2, residual stripes are visible in the NILC 
and MILCA y-maps. These stripes are mainly due to residuals 
in-scan direction systematics after subtraction of an offset for 
each Planck stable pointing period (Planck Collaboration VIII, 
2015). To study how these stripes contaminate the final y-maps 
we have decomposed them in their in-scan and cross-scan direc¬ 
tion contributions. We first convert the y-maps from Galactic to 
Equatorial coordinates for which the scan direction corresponds 
mainly to a fixed longitude value. Secondly, we apply a Galactic 
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NILC 



MILCA 




Fig. 5: Top: Standard deviation maps for the NILC (top) and 
MILCA (middle) y-maps corresponding to the inhomogenous 
noise contribution computed from the half difference of the 
half-rings maps. Bottom: Angular power spectrum of the ho¬ 
mogenous noise contribution for the NILC (orange) and MILCA 
(blue) y-maps (see main text for details). 


mask to the Equatorial y-maps and decompose them in spherical 
harmonics. Third, we select the in-scan (cross-scan) direction 
components by nullifying the spherical harmonic coefficients, 
ior i > m {i < m). Finally, we construct maps from those 
transformed coefficients and convert back to Galactic coordi¬ 
nates. We have chosen to mask the brightest 40% of the sky in 
the 857 GHz Planck map to keep Galactic ringing residuals neg¬ 
ligible and a sufficiently large fraction of the sky for the analysis. 


Figure 4 shows from left to the right the original, in-scan 
and cross-scan y-maps for NILC (top) and MILCA (bottom) in 
the southern sky region presented in Figure 3. We note that the 
y-maps look noisier as they are not filtered using the procedure 
described in the first paragraph of Sect. 6.1. The stripes are ap¬ 
parent in both the original and in-scan y-maps as expected. We 
find that the ratio of the rms of the in and cross scan maps is 
1.16 (1.17) for NILC (MILCA), consistent with residual stripe 
contamination (for Gaussian noise only maps we find a ratio of 
1). Here we have measured an overall increase of the rms of the 
maps due to stripe contamination, a more detailed estimate of 
the effect in terms of tSZ power spectrum is presented later in 
Sect. 5.2.1. 

The in-scan and cross-scan decomposition method modifies 
the cluster signal significantly due to ringing effects - see pos¬ 
itive and negative patterns around clusters in Figure 4. To ex¬ 
plore this effect we have also applied the in-scan and cross¬ 
scan decomposition method to the simulated Compton param¬ 
eter map for the detected and confirmed clusters of galaxies in 
the Planck catalogue (Planck Collaboration XXIX, 2014), which 
is presented in Section 5.3. We find that those negative and pos¬ 
itive patterns are also present in the decomposed maps and the 
ratio between the rms of the in-scan and cross-scan maps is 1 
as expected. We finally stress that the in-scan and cross-scan de¬ 
composition of the y-maps does not preserve the positiveness of 
the tSZ signal. As a consequence we can not use them to estimate 
the contamination of the stripe systematic effect in the analysis 
presented in Section 6. 

4.2. Noise distribution on they-map 

The Planck maps present a highly non-homogeneous structure 
of the noise that is mainly due to the inhomogeneous scanning 
strategy. Such complicated structure is propagated into the y- 
map and needs to be considered for further analysis. Here, we 
have chosen to describe the noise structure in the y-map by A) 
a variance map, which will capture the inhomogeneity of the 
noise, and B) the angular power spectrum of an homogenous 
noise map, which it is obtained after correcting for inhomogene¬ 
ity. The variance map and homogeneous noise power spectrum 
have been obtained using two different methods: 

1. Half-difference. The half difference of the half rings y-maps 
is used to obtain an estimate of the noise in the NILC and 
MILCA y-maps. From this half difference map we com¬ 
pute an estimate of the noise variance at lower resolution by 
squaring and downgrading it to Healpix Nside=128. The 
homogenous noise contribution is obtained by dividing the 
half difference map by the variance map after upgrading it to 
Healpix N^/je=2048. 

2. Simulations-based.lOO noise only realisation from the simu¬ 
lations described in Sect. 2.2 are used to compute an estimate 
of the variance per pixel at full resolution, which is then av¬ 
eraged to Healpix Nside=128 resolution. The homogenous 
noise contribution map is then computed as above. 

The two methods give consistent results at high Galactic lat¬ 
itudes but differ at low Galactic latitudes where we expect larger 
foreground residuals that show up in the Half-difference method. 
This is visible in the top panel of Figure 5 where we display vari¬ 
ance maps obtained from the Half-difference method for NILC 
(left) and MILCA (right). For a wide range of analyses with the 
y-map (i.e computing the radial profile, surface density or total 
fiux of a cluster, stacking of faint sources, detection of shocks, 
study of the ICM) an estimate of the overall uncertainties per 


7 


























Planck Collaboration: A map of the thermal Sunyaev-Zeldovich effect 



P5i?5oo,PSZ [arcmin^] 
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Fig. 6: Top: Comparison between the measured tSZ flux reported 
in the Planck cluster sample and that estimated directly on the 
NILC y-map for the blindly detected common sources. Bottom: 
As above but for the NILC and MILCA measured fluxes. 


pixel (including foreground contribution) is required. Thus, the 
half difference maps are released with the y-map. We also show 
in the bottom panel of Figure 5 the power spectrum of the ho¬ 
mogeneous noise contribution. We observe a signiflcant large 
angular scale (low multipoles) component in the the homoge¬ 
neous noise contribution. This low multipole component is sig- 
niflcantly larger for the NILC y-map than for the MILCA one. 

4.3. Foreground contamination and masking 

4.3.1. Galactic emission 

As discussed above, residual galactic foreground contribution is 
observed at low Galactic latitudes in both the NILC and MILCA 
y-maps, while being more important for the MILCA y-map. This 


contribution is mainly induced by thermal dust emission, as dis¬ 
cussed in details in Planck Collaboration XXIX (2014), and is 
therefore mainly associated to the Galactic plane. When dealing 
with individual objects or with the stacking of faint objects we 
recommend to account for it using the variance maps discussed 
above. For statistical analyses as those presented in Sect. 5 and 
6, speciflc masks need to be deflned and they are fully discussed 
in those sections. 

4.4. tSZ signal from resolved sources 

4.4.1. Point sources 

Point source contamination is an important issue for the cosmo¬ 
logical interpretation of the Planck Compton parameter map. 

In the reconstructed tSZ maps radio sources will appear as 
negative peaks, while infrared sources will show up as posi¬ 
tive peaks, mimicking the cluster signal. To avoid contamination 
from these sources we introduce a point source mask (PSMASK, 
hereafter). This mask is the union of the individual frequency 
point-source masks discussed in Planck Collaboration XXVIII 
(2014). The reliability of this mask was verifled by looking at 
the ID PDF of the pixel signal (as it will be discussed in Section 
6.1). We found that a more robust mask can be obtained by en¬ 
larging the mask size around the strongest radio sources in order 
to minimize their contribution. 

Alternatively, we have also performed a blind search for neg¬ 
ative sources in the y-maps using the MHW2 algorithm (Lopez- 
Caniego et al., 2006). We have detected 997 and 992 nega¬ 
tive sources for NILC and MILCA respectively. We And that 
67 and 42 (for NILC and MILCA, respectively) of those de¬ 
tected sources are not masked by the PSMASK. However, most 
of them (54 for for NILC and 36 for MILCA) are false detections 
made by the algorithm in regions surrounding very strong posi¬ 
tive sources (i.e. galaxy clusters), for those detected as additional 
negative sources, the PSMASK has been updated to account for 
them, even though the strongest are already excluded when ap¬ 
plying the 50% Galactic mask used for the cosmological analysis 
below. 

For infrared sources, estimating the efficiency of the masking 
is hampered by the tSZ signal itself. The residual contamination 
from point sources is discussed in Sects. 5.2 and 6. It is also im¬ 
portant to note that the PSMASK may also exclude some clus¬ 
ters of galaxies. This is particularly true in the case of clusters 
with strong central radio sources, such as the Perseus cluster (see 
Planck Collaboration XXIX, 2014). 

4.4.2. Blind search for clusters 

Following Planck Collaboration XXI (2014) a blind search for 
tSZ (positive) sources has been also performed on the all¬ 
sky NILC and MILCA y-maps. We use the IFCAMEX (MHW2, 
Gonzalez-Nuevo et al., 2006; Lopez-Caniego et al., 2006) and 
the single frequency matched Alter (MF, Melin et al. 2006) meth¬ 
ods. The detected sources with signal-to-noise ratio > 4 have 
been matched to the Planck cluster catalogue (PSZ2, Planck 
Collaboration XXVII, 2015). We have associated the detected 
sources to PSZ2 objects if the distance between their positions 
is smaller than 10' (the resolution of the SZ all-sky maps). The 
MHW2 algorithm flnds 1018 and 1522 candidates for NILC and 
MILCA respectively, out of which 500 and 457 correspond to 
objects present within the PSZ2 catalogue. For NILC (MILCA) 
we have 41 (25) positions that are not masked by the PSMASK 
and that do not correspond to PSZ2 clusters. However 14 (5) are 
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Fig. 7: Compton parameter MILCA (top row) and NILC (middle row) maps for a selection of PSZ2 cluster candidates with signal to 
noise ratio of 9.3, 6.2 and 4.6 from left to right. The maps are centred at the positions of the clusters in galactic coordinates, which 
are given at the bottom of the plot. Color scale is in units of y x 10“^. The bottom row presents the cluster radial profiles for MILCA 
(black) and NILC (blue). The beam profile is shown as a blue dashed line. 


located in regions excluded by the mask used to build the PSZ2 
catalogue. With the MF method we have 1472 and 1502 candi¬ 
dates for NILC and MILCA, out of which 1107 and 1096 corre¬ 
spond to objects present within the PSZ2 catalogue (867 and 835 
corresponding to already validated objects, respectively). This 
implies a good agreement between the cluster sample and the 
detected sources in the y-maps. 

This agreement is improved by taking the union of the MHW2 
and the MF catalogues. In this case, and considering only val¬ 
idated sources in the PSZ2 catalogue (1070 sources), we find 
907 and 870 matches with blind detected cluster candidates for 


NILC and MILCA, respectively. We have also performed a visual 
inspection of the NILC and MILCA y-maps at the position of 
validated PSZ2 sources for which we find no counterpart and we 
find evidence of an excess of signal. For most of these sources 
low signal-to-noise ratio and/or foreground contamination can 
explain why they are not detected blindly in the y-map. This is 
consistent with the fact that we expect blind detection methods in 
the y-map to be less sensitive than multifrequency ones (Melin 
et al., 2012). In some cases these non-detected sources are ex¬ 
tended but with a relatively high signal-to-noise (between 6 and 
12). We notice that the MHW2 and MF methods are tuned to detect 
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(a) Shapley supercluster 



- 2.0 25.0 y X 10^ 

(312.7, 30.6) Galactic 

(b) A3395-A3391 merger system 





- 2.0 25.0 y X 10^ 

(262.8, -25.2) Galactic 

(c) A339-A401 merger system 



(312.7, 30.6) Galactic 
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Fig. 8: Compton parameter maps of well known merging systems for MILCA (left column) and for NILC (right column). 
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(e) 20 < N 200 < 30 (f) 60 < N 200 < 100 



[deg] [deg] 


Fig. 10: Integrated tSZ signal as a function of cluster richness 
(left) and total mass (right). The black points correspond to 
the average signal obtained for the richness bins considered 
in Figure 9. The red line represents the corresponding best-fit 
power law (Eqs. 4 and 6 respectively). Considering z < 0.42, 
there are 13814 objects for 8 < N 200 < 10, 37250 for 10 < N 200 < 
20, 7458 for 20 < N 200 < 30, 2069 for 30 < N 200 < 40, and 1133 
for 40 < N 200 ^ 60. 


Fig. 9: 4° X 4° average maps for different ranges in N 200 from 8 
to 100. The color scale is in unit of 10“^ y. 


point-like and compact sources primarily, and thus they may fail 
for extended ones. 

We present in the top panel of Figure 6 a comparison of the 
fiux measured for these common sources in the NILC y-map 
with that derived from the PSZ2. The 75 ^ 5 ^^ correspond to the 
integrated signal within a radius equal to 5 x rsoo^. We observe 
a good agreement between the two. Most of the observed out¬ 
liers (low y-map flux with respect to the PSZ2 one) correspond 
to cluster candidates detected close to masked point sources, or 
in regions with strong Galactic contamination. We also show in 
the bottom panel a comparison of the flux measured with the 
NILC and MILCA y-maps. We find good agreement between the 
two with three outliers corresponding to sources lying either in 
a region badly affected by radio and/or infrared point sources or 
for which we observe large uncertainties in the estimated char¬ 
acteristic radius of the cluster candidate. For both maps (NILC 
and MILCA), the average signal-to-noise ratio of the recovered 
clusters with the MHW2 algorithm is ~ 10, while the signal-to- 
noise ratio of PSZ2 clusters only found with the matched filter 
technique is, as expected, lower (~ 6 ). 


^ rsoo is the cluster-centric distance at which the mean cluster density 
is equal to 500 times the critical density of the Universe 


4.4.3. Maps of selected clusters 

We have performed a more detailed analysis for some of the 
cluster candidates in the PSZ2 sample. In Figure 7 we present 
Compton parameter maps centred at the position of three of 
the newly discovered Planck clusters with signal-to-noise ratio 
of 9.3 (left), 6.2 (middle) and 4.6 (right) both for MILCA (top 
row) and NILC (middle row). In the bottom row of the figure we 
present the radial profiles of these clusters as obtained from the 
MILCA (black lines) and NILC (blue) maps. We find that the 
profiles are consistent between MILCA and NILC. For compari¬ 
son we also show as a blue dashed line the radial profile of a 1 O' 
gaussian beam. We observe that even for low-signal-to-noise 
and compact clusters we are able to detect extended emission. 
Thus, these maps could be used to extend the pressure profile 
analysis presented in Planck Collaboration Int. V (2013); Durret 
( 201 1 ) to fainter and higher redshift clusters. 

One of the major outcomes of the NILC and MILCA y- 
maps is the possibility to study diffuse faint emission between 
clusters as well as the emission in the cluster outskirts. We 
show in Figure 8 some well known merging systems includ¬ 
ing the S hap ley supercluster, and the A33 95-A33 91 and 
the A33 9-A4 01 interacting systems (see Planck Collaboration 
et al., 2013, for a detailed description). We observe that the clus¬ 
ters themselves, their outskirts and the inter cluster emission are 
well reconstructed in the NILC and MILCA y-maps, which show 
consistent results. The quality of these maps will permit analyses 
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similar to the one presented in Planck Collaboration et al. (2013) 
with a significant increase of the signal-to-noise ratio. 

4.5. SZ signal below the noise level 

We use the catalog of 132684 clusters of galaxies identified from 
the Sloan Digital Sky Survey III (Wen et al., 2012) in order to 
quantify the tSZ signal below the noise level in the y-map. This 
catalogue provides estimates of cluster redshift, richness (A^ 20 o) 
and characteristic radius (r 2 oo)- 

We focus on unresolved groups and clusters for which the 
richness is 8 < N 200 < 100. These objects are expected to be at 
signal-to-noise ratio well below 1 in the y-map and so their direct 
detection is not possible. However, their average tSZ signal can 
be detected using a stacking approach on the y-maps. Figure 9 
shows the stacked maps (obtained with the IAS stacking library, 
Bavouzet 2008; Bethermin et al. 2010) for 6 richness intervals 
with N 200 going from 8 to 100, on patches of 4° x 4°, for the full 
sky NILC y-map. For the stacking we exclude all the positions 
for which point sources are found within a radius of 10 arcmin 
from the center of the object. The noise of the stacked maps 
scales as expected with 1 / yfN and the average signal increases 
for increasing richness. Then, when considering a sufficiently 
large number of objects, we are able to significantly detect the 
average SZ signal, even for small groups (Fig. 9). We obtain 
consistent results for the MILCA and NILC y-maps. 

For the N 200 intervals reported in Fig. 9 we have estimated 
the total stacked fluxes (^ 5 ^ 5 ^^) as the integrated signal within a 
radius 5 x rsoo, with r^m = 0 . 7 (r 2 oo), <^ 200 > being the mean of 
the r 2 oo reported by Wen et al. (2012) for clusters belonging to 
each considered subsample. The average values and associated 
errors have been obtained with a bootstrap approach. For this, we 
have constructed and stacked cluster samples obtained by ran¬ 
domly replacing sources from the original sample, so that each 
of them contains a number of clusters equal to the initial one. 
By repeating the process several times, we have determined the 
statistical properties of the population being stacked. In Figure 
10 (left panel) we report 75^500 as a function of richness (A^ 20 o) 
for all the objects with z < 0.42. In fact higher redshift clusters 
present in this catalogue may have a biased lower richness be¬ 
cause of incompleteness of member galaxies, as detailed in Wen 
et al. (2012). Following Planck Collaboration et al. (2011), we 
have fitted a power law of the form 



with E{z)^ = and Da the angular diameter distance 

of the cluster. 

By considering only A ^200 ^ 60 (for which the number 
of clusters in each richness bin is > 1000 ), we find a slope 
a = 2.21 + 0 . 10 , which is consistent with the scaling obtained by 
Planck Collaboration et al. (2011). We have verified that the scal¬ 
ing obtained is insensitive to wether we limit our cluster sam¬ 
ple to high Galactic latitude objects or not, and that a change in 
the integration radius only affects the normalization 720 - With 
= 75 , 3 ,,/1.79 (Planck Collaboration XXIX, 2014), we ob¬ 
tain 720 = (8.07 + 0.41) X 10“^arcmin^. However, for indepen¬ 
dent cluster datasets, different choices are made for the fainter 
magnitude of the member galaxies contributing to A^ 200 - Then 
this may affect the richness associated to an object of given mass 
and complicate the comparison of the constant term 72 o for anal¬ 
ysis that are based on different cluster datasets. 


The mass-richness scaling relation, 

M2oo = Mo(^)^ (5) 

can then be used to also derive the 7500 -Mtot scaling (see 
Figure 10, right panel), by fitting a power low of the form 

Ford et al. (2014) compares different results for the scaling of 
the mass with richness, in the form of Eq. 5. In particular they 
discuss those obtained by Wen et al. (2012), for a subsample 
of clusters with already known masses, and by Covone et al. 
(2014), from weak lensing mass measurements of 1,176 clus¬ 
ters of the catalogue used here. Assuming p = 1.2 + 0.1 and 
Mo = (1.1 ± 1.1) X 10^"^ Mq (Wen et al., 2012; Covone et al., 
2014, and a ten percent error) in Eq. 5, we compute via a 
Monte Carlo approach the average M 200 and associated uncer¬ 
tainties for each bin. From this we find B = 1.92 + 0.42 
for the 75^3 qq-M2oo scaling (Eq. 6). This is consistent with the 
self-similar value (5/3) and with what is expected for this mass 
range according to simulations (Sembolini et al., 2014). Through 
a weak lensing analysis of > 18000 clusters in the CFHTLenS 
(Canada-France-Hawaii Telescope Lensing Survey), Ford et al. 
(2014) obtain a steeper calibration for the mass-richness rela¬ 
tion, with p = 1.4 + 0.1. For an higher p, the slope B 

gets lower {B = 1.66 + 0.37), but still consistent with the value 
obtained before. Again comparison of Mq and 7o is complicated 
by the different choices of the fainter limit on galaxies contribut¬ 
ing to A/ 200 • By exploring the range 0.2 < z < 0.9, Ford et al. 
(2014) finds no evolution with redshift. Thus, we do not expect 
our result to be biased by the z < 0.42 cut. 

5. Angular power spectrum 

We now consider the validation of the Compton parameter maps 
at the power spectrum level. 

5.1. Power spectrum estimation 

To estimate the power spectrum of the tSZ signal we use the 
XSPECT method (Tristram et al., 2005) initially developed for 
the cross-correlation of independent detector maps. XSPECT 
uses standard MASTER-like techniques (Hivon et al., 2002) to 
correct for the beam convolution and the pixelization, as well 
as the mode-coupling induced by masking foreground contami¬ 
nated sky regions. 

In the following, all the spectra will use a common multi¬ 
pole binning scheme, which was defined in order to minimize 
the correlation between adjacent bins at low multipoles and to 
increase the signal-to-noise at high multipole values. We will 
also only consider cross angular power spectra between the y- 
maps reconstructed from the first (F) and second (L) halves of 
the data. This allows us to avoid the bias induced by the noise in 
the auto angular power spectrum, the correction of which would 
require a large number of noise simulations. The cross angu¬ 
lar spectra are named in the following NILC F/L and MILCA 
F/L. Error bars in the spectrum are computed analytically from 
the auto-power and cross-power spectra of the pairs of maps, 
as described in Tristram et al. (2005). We do not consider here 
higher order terms to the power spectrum variance. All of our 
Compton parameter maps assume a circular Gaussian beam of 
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Fig. 11: Angular cross-power spectra of the Planck NILC F/L (left) and MILCA F/L (right) reconstructed Compton parameter maps 
for different Galactic masks corresponding to 30 % (cyan), 40 % (green), 50 % (blue) and 60 % (pink) of the sky. For comparison 
we also show MILCA-NILC F/L (red) and NILC7 0 VfL (black) on 50 % of the sky. See text for details. 


10' FWHM. The additional filtering at large angular scales in 
the MILCA Compton parameter maps is also accounted for and 
deconvolved. 

5.2. Foreground contamination 

In Planck Collaboration XXI (2014) we have identified the dom¬ 
inant foreground contributions to the angular power spectrum of 
the reconstructed y-maps using the FFP6 simulated data set. We 
have repeated here the same analysis on updated simulated data 
sets (see Sect. 2.2) for which the description of foreground com¬ 
ponents has been improved according to our current knowledge 
(Planck Collaboration XII, 2015). We find no major difference 
between the two analyses and therefore we do not repeat the dis¬ 
cussion here. At large angular scales the dominant foreground 
contribution is the Galactic thermal dust emission as we have 
already discussed in Section 4.3.1. At intermediate and small 
angular scales the major contribution comes from the clustered 
CIB emission. Equally important at small angular scales are the 
residual contribution from radio and IR sources. 

5.2.1. Low-multipole contribution 

Assuming that at large angular scales the Compton parameter 
maps are mainly affected by diffuse Galactic dust emission, we 
have tested several Galactic masks by imposing fiux cuts on the 
Planck 857 GHz channel intensity map. In particular we inves¬ 
tigated masking out 30%, 40%, 50% and 60% of the sky. The 
edges of these masks have been apodized to limit ringing effects 
on the reconstruction of the angular power spectrum. Figure 11 
shows the power spectra for NILC (left) and MILCA (right). 
We observe that the MILCA F/L large angular scales power de¬ 
creases when imposing a more severe masking (larger fraction 
of the sky is masked out). A similar effect, but significantly 
smaller, is also observed for the MILCA F/L cross power spec¬ 
trum. The MILCA y-map is significantly more contaminated by 
thermal dust emission at multipoles below 30. However, we also 
observe for NILC F/L an extra noise contribution at large angu¬ 
lar scales as discussed in Section 4.2. These two problems limit 
the reliability of the results at multipoles below 30. 

To extend the measurement of the angular power spectrum 
of the tSZ emission to multipoles below 30 we have considered 
two options: 1) as in Planck Collaboration XXI (2014) we ap¬ 


ply a more severe galactic mask (30 % of the sky is masked) 
before the computation of the NILC weights to produce the y- 
map (NILC7 0), and 2) compute the cross-correlation of NILC 
and MILCA y-maps. Considering 50% of the sky we show in 
Figure 11 the angular power spectrum of the NILC7 0 y-map 
(black) and the cross spectrum of the NILC first half and MILCA 
second half ( NILC-MILCA F/L) y-maps. We observe that the 
two are compatible within error bars in the multipole range 
10 < I < 1500. As expected the cross spectrum NILC-MILCA 
F/L shows larger error bars. 

Although the NILC7 0 y-map seems to be the best choice 
in terms of power spectrum estimation, it results in a significant 
reduction of the available sky area for other kind of studies as 
those presented in Section 4. Furthermore, it is difficult to ac¬ 
curately estimate the ultimate residual foreground contribution 
at very large angular scales. Because of this and to preserve the 
coherence of the delivered products and the analysis presented 
in this paper, we have chosen the cross angular power spectrum 
of NILC-MILCA F/L as a baseline. This is obviously a more 
conservative choice in terms of noise induced uncertainties. The 
NILC-MILCA F/L cross angular power spectrum bandpowers 
and uncertainties are further discussed in Sect. 7. Using the in¬ 
scan and cross-scan y-maps presented in Section 4.1 we find that 
stripe contamination accounts for less than 10% of the total sig¬ 
nal in the NILC-MILCA F/L cross angular power spectrum. We 
have enlarged the error bars to account for this systematic effect. 

5.2.2. High-multipole contribution 

At small angular scales the measured tSZ power spectrum is af¬ 
fected by residual foreground contamination coming from clus¬ 
tered CIB emission as well as radio and IR point sources. They 
show up in the MILCA-NILC F/L cross power spectrum (see 
Figure 1 1) as an excess of power at large multipoles. 

To deal with those we adopt the same strategy as in Planck 
Collaboration XXI (2014). We define physically motivated mod¬ 
els of the angular power spectrum of the foreground components 
for each observation channel, including cross correlations be¬ 
tween channels. In contrast to Planck Collaboration XXI (2014) 
we also account for the cross correlation between the clustered 
CIB and the tSZ emission. A detailed description of this cross 
correlation as well as of the clustered CIB model is presented in 
Planck Collaboration XXIII (2015). For the radio and IR point 
source models we refer to Planck Collaboration XXI (2014). 
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Fig. 12: Comparison of the tSZ angular power spectrum estimated from the cross-power-spectrum of the NILC (left) and MILCA 
(right) F/L maps (black) with the expected angular power spectrum of the confirmed clusters in the Planck Cluster Sample (green 
line). The angular cross-power spectrum between the NILC A and MILCA Compton parameter maps and the simulated detected 
cluster map is shown in red. Tthe correlation between the reconstructed y-map and the simulated detected cluster map, to which an 
arbitrary rotation has been applied, is plotted in grey. 


Using the models described above we compute Gaussian re¬ 
alizations of the foreground contribution for each HFI frequency 
channel between 100 and 857 GHz. Note that the LFI channels 
are only used at large angular scales. We apply the MILCA or 
NILC weights to these simulated maps. From these simulations 
we find that the cross correlation between the CIB and tSZ con¬ 
tribution can be neglected to first order with respect to the others 
and it will therefore not be considered hereafter. Uncertainties 
on the parameters describing the foreground models were also 
propagated using simulations. We find that the clustered CIB 
model uncertainties might be as large as 50% in amplitude. In 
addition, we notice that the amplitude of the point source mod¬ 
els can vary significantly with the point source mask applied. 
These uncertainties are taken into account hereafter. The ampli¬ 
tude of the residual foreground models are jointly fitted with the 
cosmology-dependent tSZ model, as detailed in Sect. 7.1. 

5.3. Contribution of resolved clusters to the tSZ power 
spectrum 

We simulate the expected Compton parameter map for the de¬ 
tected and confirmed clusters of galaxies in the Planck cata¬ 
logue (Planck Collaboration XXIX, 2014) from their measured 
integrated Compton parameter, 751 - 5 ^^ and redshift, z. We assume 
hydrostatic equilibrium and an Arnaud et al. (2010) pressure pro¬ 
file. The green solid line in Fig. 12 shows the power spectrum of 
this simulated map. Figure 12 also shows the cross-power spec¬ 
tra of the NILC and MILCA F/L maps (in black). We also com¬ 
pute the cross-power spectrum of the simulated cluster map and 
the Planck reconstructed NILC and MILCA Compton parameter 
maps. This is shown in red in the figure. Here again, the signal 
is consistent with the expected power spectrum of the confirmed 
Planck clusters of galaxies. These results show that the tSZ sig¬ 
nal from clusters is preserved in the y-maps. 

6. Higher order statistics 

The power spectrum analysis presented above only provides in¬ 
formation on the 2-point statistics of the Compton parameter dis¬ 
tribution over the sky. An extended characterization of the field 
can be performed by studying the higher-order moments in the 


ID PDF of the map, or by measuring 3-point statistics, i.e., the 
bispectrum. 

6.1. ID PDF analysis 

Following Planck Collaboration XXI (2014) we perform an 
analysis of the ID PDF of the NILC and MILCA reconstructed 
Compton parameter maps. For the tSZ effect we expect an asym¬ 
metric distribution with a significantly positive tail (Rubino- 
Martm & Sunyaev, 2003). We thus focus on the asymmetry of 
the distribution and its unnormalized skewness. First, we filter 
the maps in order to enhance the tSZ signal with respect to 
foreground contamination and noise. We apodise the combined 
PSMASK and Galactic mask to avoid residual point source ring¬ 
ing. We follow the approach of Wilson et al. (2012) and use 
a filter in harmonic space. For each multipole the filter is 
constructed as the ratio between the angular power spectrum of 
the expected tSZ signal obtained from the simulations in 

Sect. 2.2) and the power spectrum of the half-difference y maps 
(C^, see Sect. 4.2) such that = Cf^lC^. We smooth this 
filter function in multipole space using a 21 -point square ker¬ 
nel and normalize it to one. Notice that this filter only selects 
the multipole range for which the tSZ signal is large with re¬ 
spect to the noise, and thus, it does not modify significantly the 
non-Gaussianity properties of the y-maps. Furthermore, we have 
found that the filter used here behaves better than the more tradi¬ 
tionally used Wiener filter, as it is less affected by point-source 
ringing. Following this procedure, the ID PDF of the filtered 
Compton parameter maps, P(y), is computed from the histogram 
of the pixels. As for the power spectrum, several Galactic masks 
have been considered in order to tests the robustness of the re¬ 
sults. 

Figure 13 shows the ID PDF for the NILC (left) and MILCA 
(right) Compton parameter map in red when masking 50% of 
the sky. Each of them corresponds to the convolution of the ID 
PDF of the different components in the map: the tSZ effect; fore¬ 
grounds; and the noise. The ID PDF of the NILC and MILCA 
y-maps clearly show three distinct contributions: a Gaussian cen¬ 
tral part that is slightly wider than the noise contribution, as ex¬ 
pected from the half-difference map ID PDF (black curve); a 
small negative tail, corresponding most likely to residual radio 
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Fig. 13: ID PDF of the Planck y-map before (red) and after (orange) masking the PSZ2 clusters, and of the half-difference map 
(black) for the NILC (left) and MILCA (right) methods. We also show for comparison the ID PDF of the simulated PSZ2 cluster 
map (dark blue). 


sources; and a positive tail corresponding mainly to the tSZ sig¬ 
nal as observed for the PSZ2 cluster simulated map (dark blue). 
In comparison to the Planck Collaboration XXI (2014) we find 
now a better agreement between the NILC and MILCA results. 
This is mainly due to the improved masking of radio sources pre¬ 
sented in Section 4.4.1. Finally we also show the ID PDF of the 
reconstructed y-maps after masking the PSZ2 clusters (orange). 
We observe that most of the tSZ is removed, indicating that the 
ID PDF of the y-map is dominated by detected clusters. 

A simple analysis of the measured ID PDF can be performed 
by considering the asymmetry of the distribution: 

r+oo 

P(y)dy- P(y)dy, (7) 

*Jyp *J—oo 

where yp is the peak value of the normalized PDF P(y). In ad¬ 
dition, the non-Gaussianity of the positive tail can be quantified 
by 

^+oo 

A= [P(y)-G(y)]dy, (8) 

dyp 

with G(y) the expected PDF if fluctuations were only due to 
noise. The latter can be obtained from the half-difference y- 
maps. 

Masking 60 % of the sky, we find A = 0.218 and A = 0.10 for 
the NILC Compton parameter map. Equivalently, for the MILCA 
Compton parameter map we find A = 0.223 and A = 0.11. These 
results are in agreement and consistent with a positive tail in the 
ID PDF, confirming the tSZ nature of the signal. The agreement 
between the NILC and MILCA results degrade when reducing 
the masked area as a consequence of a larger foreground contri¬ 
bution in the MILCA y-map. See Hill et al. (2014) for a similar 
analysis conducted on ACT data. 

6.2. Bispectrum 

Since the SZ signal is non-Gaussian, significant statistical 
information is contained in the bispectrum, complementary 
to the power spectrum (Rubino-Martm & Sunyaev, 2003; 
Bhattacharya et al., 2012). Results on SPT data have been ob¬ 
tained by Crawford et al. (2014) and George et al. (2014). The 
bispectrum also provides an additional statistic to assess the 
compatibility of the NILC and MILCA reconstructed Compton 
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Fig. 14: Bispectra of the NILC (green) and MILCA (red) y-maps 
for four different configurations (equilateral, orthogonal, fiat and 
squeezed), compared with the bispectrum of the projected map 
of the PSZ2 clusters (blue). ±\(t uncertainties are indicated as 
black dotted lines. 


parameter maps, as well as their reliability in terms of fore¬ 
ground contamination. 

We therefore analyze the bispectra of the NILC and MILCA 
maps. The estimation method is essentially the same as for 
the 2013 results (Planck Collaboration XX, 2014), briefiy re¬ 
capped here. We use the binned bispectrum estimator described 
in Bucher et al. (2010) and Lacasa et al. (2012), which is also 
used for the Planck primordial non-Gaussianity analysis (Planck 
Collaboration XXIV, 2014; Planck Collaboration XVII, 2015). 
We use a multipole bin size AT’ = 64 and a maximum mul¬ 
tipole Tmax = 2048 for the analysis, working at a resolution 
A^side = 1024 to reduce computing time. The NILC and MILCA 
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Fig. 15: NILC - MILCA F/L cross-power spectrum (black) com¬ 
pared to the power spectra of the physically motivated fore¬ 
ground models. The considered foregrounds are: clustered CIB 
(green line); infrared sources (cyan line); and radio sources (blue 
line). Additionally, the best-fit tSZ power spectrum model pre¬ 
sented in Sect. 7.1 is also plotted as a solid red line. 


maps are masked with the combination of PSMASK, described 
in Sect. 4.4.1, and a Galactic mask at 50, 60 or 70%, described 
in Sect. 5.2.1 (in the rest of this section we will simply denote as 
X% mask the combination of PSMASK and the Galactic mask 
at X%). The best-fit monopole and dipole outside the mask are 
finally removed before estimation. 

An important part of the pipeline is then to correct for the 
bias introduced by masking. To this end, we compute the ratio of 
the full-sky and masked sky bispectra, on highly non-Gaussian 
simulations with a tSZ-like bispectrum and a 10' resolution. This 
ratio is used to correct the measured bispectra and fiag unre¬ 
liable (^ 1 ,^ 2 , ^ 3 ) configurations. Specifically we fiag configura¬ 
tions where the ratio is different by > 25% from the naive ex¬ 
pectation /sky B{ix)B{i 2 ) ^(^ 3 ), where B{i) is the Gaussian 10' 
beam. 

For both NILC and MILCA, we find that the bispectra com¬ 
puted on the 50, 60 and 70% masks are consistent. This indi¬ 
cates that there is no detectable residual galactic contamination 
in these bispectra. However we did find Galactic contamination 
on less aggressive Galactic masks, specifically positive Galactic 
dust. As Galactic dust is highly non-Gaussian, we warn the use 
against the measurement of higher order statistics using Galactic 
masks smaller than 50%. In the following we adopt the 50% 
mask as baseline, as it leaves the most sky available for estima¬ 
tion and minimizes masking effects in the measurement. 

Figure 14 shows the obtained bispectra as a function of mul¬ 
tipole for the NILC (green) and MILCA (red) Compton param¬ 
eter maps. We observe a good agreement between the bispectra 
of the two maps, and the bispectral behaviour is consistent with 
that expected from a tSZ signal (see e.g. Lacasa, 2014, chap¬ 
ter 5). We furthermore compare these measurements with the 
bispectrum of the simulated map for the PSZ2 clusters, which 
is presented in blue in Fig. 14. We observe a good agreement 
between the bispectra of the NILC and MILCA and that of the 
PSZ2 clusters. We therefore conclude that the observed bispec¬ 
trum in the y-map is dominated by detected clusters. 

Finally, in Fig. 14 are shown the ±\(t uncertainties of the 
measurements, in black dotted lines. The error bars were com¬ 
puted in a similar manner to that of the 2013 results (Planck 


Table 2: Marginalized bandpowers of the angular power spec¬ 
trum of the Planck tSZ Compton parameter map (in dimension¬ 
less (AT/r)^ units), statistical and foreground errors, and best-fit 
tSZ power spectrum and number counts models (also dimension¬ 
less). 


f • 

^min 

^max 

4ff 

£(£ + l)Q/2;r 

[lO'V] 

^stat 

[lo'V] 

CTfg 

[io‘V] 

Best-fit 

[io‘V] 

9 

12 

10.0 

0.00506 

0.00629 

0.00002 

0.00726 

12 

16 

13.5 

0.00876 

0.00615 

0.00007 

0.00984 

16 

21 

18.0 

0.01353 

0.00579 

0.00015 

0.01320 

21 

27 

23.5 

0.02946 

0.00805 

0.00021 

0.01737 

27 

35 

30.5 

0.02191 

0.00522 

0.00053 

0.02274 

35 

46 

40.0 

0.02744 

0.00464 

0.00109 

0.03008 

46 

60 

52.5 

0.04093 

0.00468 

0.00172 

0.03981 

60 

78 

68.5 

0.04227 

0.00429 

0.00320 

0.05236 

78 

102 

89.5 

0.06463 

0.00454 

0.00567 

0.06901 

102 

133 

117.0 

0.10738 

0.00562 

0.00969 

0.09102 

133 

173 

152.5 

0.12858 

0.00594 

0.01889 

0.11956 

173 

224 

198.0 

0.15696 

0.00611 

0.02895 

0.15598 

224 

292 

257.5 

0.21738 

0.00687 

0.04879 

0.20306 

292 

380 

335.5 

0.28652 

0.00824 

0.08374 

0.26347 

380 

494 

436.5 

0.36682 

0.00958 

0.13524 

0.33848 

494 

642 

567.5 

0.42666 

0.01242 

0.19500 

0.42930 

642 

835 

738.0 

0.53891 

0.01645 

0.27718 

0.53577 

835 

1085 

959.5 

0.71103 

0.02402 

0.37576 

0.65454 

1085 

1411 

1247.5 

0.82294 

0.04172 

0.55162 

0.77885 


Collaboration XX, 2014), see Appendix A.3 for a more detailed 
discussion. 

With a detection per configuration at an average significance 
of 3.5(T, and a total significance of ~ 60(T, the Planck data thus 
provide a high quality measurement of the non-Gaussianity of 
the thermal Sunyaev-Zel’dovich signal, with undetectable con¬ 
tamination from foregrounds. 

7. Cluster physics and cosmology 

7 . 1. Power spectrum analysis 

7.1.1. tSZ power spectrum modelling 

As a measure of structure growth, the tSZ power spectrum can 
provide independent constraints on cosmological parameters. 
As shown by Komatsu & Seljak (2002), the power spectrum 
of the tSZ effect is highly sensitive to the normalization of the 
matter power spectrum, commonly parameterized by the rms of 
the z = 0 mass distribution on 8 Mpc scales, erg, and to the 
total amount of matter Qm- We expect the tSZ power spectrum 
to also be sensitive to other cosmological parameters, e.g., the 
baryon density parameter Qb, the Hubble contant Hq, and the 
primordial special index n^. For reasonable external priors on 
those parameters, however, the variations are expected to be 
negligible with respect to those introduced by changes in Qm 
and (Tg and so they are not considered here. 

Following Planck Collaboration XXI (2014) we consider 
here a two-halo model for the tSZ power spectrum, which is 
fully described in Appendix A. 1 . This model accounts for both 
intra-halo (1-halo term) and inter-halo (2-halos) correlations. 
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Fig. 16: 2D and ID likelihood distributions for the combination of cosmological parameters (T8(Dni/0.28)^/^, and for the foreground 
parameters ARad.ps. ^cm and Ajr.ps. We show the 68.3% and 95.4% C.L. contours. The red and black contours correspond to a fixed 
mass bias of 0.2 and 0.4 respectively. 


Following Eq. (A.6), the tSZ spectrum is computed using the 
2-halo model, the Tinker et al. (2008) mass function, and the 
Arnaud et al. (2010) pressure profile. In particular, we use the 
numerical implementation presented in Taburet et al. (2009, 
2010, 201 1), and integrating in redshift from 0 to 3 and in mass 
from 10^^ Mq to 5 x 10^^ M©. Our model allows us to com¬ 
pute the tSZ power spectrum at the largest angular scales. It 
is consistent with the tSZ spectrum presented in Efstathiou & 
Migliaccio 2012 (EM 12), which was used as a template in the 
CMB cosmological analysis in Planck Collaboration XVI (2014) 
and Planck Collaboration XI (2015). We also include the mass 
bias parameter, b, which accounts for bias between the obser- 
vationally deduced and true (M^^“^) mass of the clus¬ 

ter (see Planck Collaboration XX, 2014, for details) such that 

Mobs ^ 

Cosmological parameter results are very sensitive to the 
mass bias and in particular we expect erg and Dm to be strongly 
degenerated with b (Planck Collaboration XXI, 2014). By con¬ 
trast to Planck Collaboration XXIV (2015), for which external 
priors in the mass bias have been used, we consider here only 
two distinct values: b = 0.2 and b = 0.4. The former corresponds 
to the average value that numerical simulations seem to indicate 


(Pig. A.2 in Planck Collaboration XX, 2014). The latter value 
for the bias alleviates the inconsistency with the constraints de¬ 
rived from the analysis based on primary CMB anisotropies (see 
Planck Collaboration XX, 2014). 

This value is however larger than that obtained by mass com¬ 
parison of clusters present in both the Planck cosmology sample 
(Planck Collaboration XX, 2014) and the Weighing the Giants 
(WtG von der Linden et al., 2014) project. Even if studies based 
on lensing mass measurements still provide different and incon¬ 
sistent results for the cluster mass calibration, their number and 
their accuracy has incredibly improved in the very recent past 
{i.e. Mahdavi et al., 2013 for the CCCP project, Umetsu et al., 
2014 for CLASH, Israel et al., 2014 for 400d WL, Lord et al. 
(2014) for CFHTLenS, Gruen et al., 2014 for WISCy) and they 
are expected to provide useful information. 

7.1.2. Maximum likelihood analysis 

As in Planck Collaboration XXI (2014), cosmological con¬ 
straints are obtained from a fit of the tSZ power spectrum. As 
discussed in Section 5.2 we take the NILC-MILCA P/L cross¬ 
power spectrum (black dots in Pigure 15) as a reference and limit 
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the analysis to 50% of the sky to minimize foreground residuals. 
In terms of astrophysical signal we consider a four-component 
model: tSZ, clustered CIB, radio point sources, and infrared 
point sources. We restrict the analysis to multipoles i > 10 so 
that we can neglect the residual thermal dust contamination (see 
Section 5.2.1). For i > 2000 the total signal in the tSZ map 
is dominated by correlated noise which is also accounted for in 
the fit. Because of this correlated noise and the expected high 
value of foreground contamination we limit the fit to multipoles 
i < 1411. Finally, the observed y-map power spectrum, C^, is 
modelled as: 

Cf = Cf^(Qm, trg) + Acib C™ + (9) 

Air + ARad + Acn 

Here cts) is the tSZ power spectrum (in red in Fig. 

15), C™ is the clustered CIB power spectrum (in green), and 
and are the infrared and radio source power spectra, 
respectively (in cyan and in blue). is an empirical model for 
the high multipole correlated noise. 

Foreground contamination is modelled following Sect. 5.2.2. 
As discussed there, the main uncertainties in the residual power 
spectrum translate into up to 50% uncertainty in the amplitude 
of the clustered CIB. We have not considered in the analysis 
the CIB-tSZ cross-correlation that was proved to be negligible 
in terms of power spectrum. The amplitude of the IR and ra¬ 
dio point-source contribution will depend very much on the ex¬ 
act Galactic mask used for the analysis. However, we expect the 
shape of the their power spectra to remain the same. We thus al¬ 
low for a variation of the normalization amplitudes for the clus¬ 
tered CIB, Acib, and for the point sources. Air and A^ad- 

We assume a Gaussian approximation for the likelihood 
function. Best-fit values and uncertainties are obtained using an 
adapted version of the CosmoMC algorithm (Lewis & Bridle, 
2002). Only erg and Qm are allowed to vary. All other cosmo¬ 
logical parameters are fixed to their best-fit values as obtained 
in Table 2 of Planck Collaboration XVI (2014). The normal¬ 
ization amplitudes, Acib, Aj-ad and Air, considered as nuisance 
parameters, are allowed to vary between 0 and 3. For the range 
of multipoles considered here, the tSZ angular power spectrum 
varies like oc The results are thus presented in terms 

of this parameter combination. 

7.2. Best-fit parameters and tSZ power spectrum 

Figure 16 presents the 2D and ID likelihood distributions for 
the cosmological parameter combination erg(Dm/0.28)^^^ and 
for the foreground nuisance parameters. We present the results 
obtained assuming a mass bias of 0.2 (black) and 0.4 (red). We 
obtain very similar values for the nuisance parameters in both 
cases. In particular the best-fit values for a mass bias of 0.2 are 
Acib = = 0.01:°;™ and A,r = 1.97:0;20. However, 

there is a significant shift in the value of crg(Qn,/0.28)^^^ as one 
would expect (Planck Collaboration XX, 2014). In the case of 
a mass bias of 0.2 we have (Tg(Qni/0.28)^/^ = 0 . 80 ^qq 3 , while 
for a mass bias of 0.4 we have (Tg(Qni/0-28)^/^ = 0.90 ^qq 3. 
Notice that these values are obtained in a specific framework: 
all other cosmological parameters being fixed and a fiducial 
fixed model used for the signals. Relaxing this framework would 
likely weaken the constraints presented here as discussed below. 

Figure 15 shows the NILC-MILCA F/L angular cross-power 
spectrum before correcting (black dots) for foreground contri¬ 
bution. We also show the best-fit foreground models: clustered 
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Fig. 17: NILC - MILCA F/L cross-power spectrum after fore¬ 
ground subtraction (red points), compared to the Atacama 
Cosmology Telescope (ACT; cyan dot) and the South Pole 
Telescope (SPT; orange, George et al., 2014) power spectrum 
estimates. The black line shows the tSZ power spectrum tem¬ 
plate (EM 12 Efstathiou & Migliaccio, 2012) used in the Planck 
CMB cosmological analysis (Planck Collaboration XVI, 2014; 
Planck Collaboration XI, 2015) with its best fit amplitude Atsz 
(Planck Collaboration XI, 2015), the grey region allows compar¬ 
ison with the 2 (T interval. 


CIB (green line), and radio (blue line) and IR (cyan line) point 
sources. The statistical (thick line) and total (statistical plus fore¬ 
ground, thin line) are also shown. The best-fit tSZ power spec¬ 
trum is presented as a solid red line. We conclude that the NILC- 
MILCA F/L angular cross-power spectrum is dominated by tSZ 
for multipoles ^ < 700, and by foreground contribution for mul¬ 
tipoles i > 1200. We also note that for the best-fit model the 
radio point-sources contribution seems to be negligible with re¬ 
spect to the IR one. This is not a physical result and it is most 
probably explained by the strong degeneracy observed between 
the radio and IR point-source amplitude (see Fig. 16). 

Finally we present in Figure 17 the NILC-MILCA F/L angu¬ 
lar cross-power spectrum after correcting for foreground contri¬ 
bution (red dots). Uncertainties account for statistical and sys¬ 
tematic errors as well as for uncertainties in the foreground 
subtraction. The marginalized bandpowers and uncertainties are 
also presented in Table 2. We note that foreground induced un¬ 
certainties dominate at multipoles i > 100. Bandpowers for 
the best-fit model for the angular tSZ power spectrum are also 
given for comparison. We also show in the Figure 17 tSZ power 
spectrum estimates at high multipoles obtained in CMB oriented 
analyses by the Atacama Cosmology Telescope (ACT; cyan dot) 
and the South Pole Telescope (SPT; orange, George et al., 2014). 
The black line shows the tSZ power spectrum template (EM 12 
Efstathiou & Migliaccio, 2012) used in the Planck CMB cos¬ 
mological analysis (Planck Collaboration XVI, 2014; Planck 
Collaboration XI, 2015) assuming the best-fit amplitude Atsz 
and the grey region 2cr uncertainties from Planck Collaboration 
XI (2015). We observe that the amplitude of the tSZ signal found 
in this paper is consistent with the high multipole based measure¬ 
ments. 
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Fig. 18: tSZ power spectrum for existing models in the litera¬ 
ture. NILC-MILCA F/L cross-power spectrum after foreground 
correction (black dots) compared to the Atacama Cosmology 
Telescope (ACT; cyan dot) and the South Pole Telescope (SPT; 
orange, George et al., 2014) power spectrum estimates. We 
also show the tSZ power spectrum models from hydrodynamic 
simulations (Battaglia et al., 2012, blue), from A^-body simu¬ 
lations plus semi-analytical dust gas models (Trac et al., 2011, 
cyan; TB02), and from analytical calculations (Shaw et al. 2010, 
green). 


7.2.1. Cluster physics dependence 

As discussed in Planck Collaboration XXI (2014), we also ex¬ 
pect the tSZ power spectrum amplitude to be sensitive to the 
physics of clusters of galaxies. To explore this dependence we 
have considered a set of predicted tSZ spectra for various phys¬ 
ical models. In Fig. 18 we compare these models to the fore¬ 
ground cleaned Planck tSZ power spectrum derived above (grey 
dots), as well as to the Atacama Cosmology Telescope (ACT; 
cyan dot) and the South Pole Telescope (SPT; orange, George 
et al., 2014) power spectrum estimates. We consider the predic¬ 
tions derived from hydrodynamical simulations (Battaglia et al., 
2010; Battaglia et al., 2012, blue), from Wbody simulations plus 
semi-analytical models (Trac et al., 2011, cyan; TB02) and from 
analytical calculations (Shaw et al. 2010, green). These models 
were originally computed for the set of cosmological parameters 
in Hinshaw et al. (2012) with erg = 0.8 and have been rescaled 
in amplitude to our best-fit value for We note that there 

is some dispersion in the predicted amplitudes and shapes of 
the tSZ power spectrum. These differences refiect the range of 
methodologies and assumptions used both in the physical prop¬ 
erties of clusters and in the technical details of the computation. 
The latter includes differences in the redshift ranges and also in 
the mass intervals probed by the limited sizes of the simulation 
boxes of the hydrodynamical simulations. Analytical predictions 
are also sensitive to the model ingredients, such as the mass func¬ 
tion, mass bias and scaling relations adopted. 

We see from Fig. 18 that the models presented above (the tSZ 
template for CMB analyses, plus the Battaglia et al. 2012, Shaw 
et al. 2010 and TB02 models) provide reasonable fits to the data 
for multipoles above 200. For lower multipoles the Shaw et al. 
2010 and TB02 models are not consistent with the data. 

We have also performed a simplified likelihood analysis to 
evaluate the uncertainties in cosmological parameters induced 
by the uncertainties in the modelling of the cluster physics. We 


replace our own model of the tSZ power spectrum by the models 
discussed above and recompute (T8(^lm/0.28)^/^, Acm, ARad and 
Air from a simple linear fit to the NILC-MILCA F/L cross-power 
spectrum. In the case of mass bias of 0.2, we obtain values for 
(T8(^lm/0.28)^/^ between 0.77 and 0.80, which lie within the 1 cr 
uncertainties (0.03) presented above. 

In the case of our fiducial model (see Appendix A.l) we can 
also consider uncertainties in the parameters describing the scal¬ 
ing relations allowing us to relate the observed tSZ fiux to the 
mass of the cluster for a given redshift. Following Eq. (7) in 
Planck Collaboration XXVIII (2014) the main parameters to be 
considered are the mass bias b, the overall amplitude F* and the 
scaling slope (3. As discussed above the mass bias is fully degen¬ 
erate with (Tg. Similar conclusions can be drawn for F*, which is 
expected to be known at the percent level (see Table 1 in Planck 
Collaboration XXVIII (2014)) and therefore it is subdominant 
with respect to the uncertainties in the mass bias. Although the 
uncertainties in the slope of the scaling relation are relatively 
large, we have checked that they lead to negligible uncertainties 
on cosmological parameters. 

7.3. Higher order statistics 
7.3.1. Skewness measurements 

The skewness of the ID PDF distribution, 
Jy^P(y)dy/{fy^P(y)dy) can also be used to derive 
constraints on cosmological parameters. Following Wilson 
et al. (2012); Planck Collaboration XXI (2014) we have chosen 
a hybrid approach, by computing the skewness of the filtered 
Compton parameter maps outside the 50% sky mask. In 
particular, we have computed the skewness of the Planck data 
Compton parameter maps (y^), and of the half-difference maps 

Using the models presented in Sect. A we can show that the 
unnormalized skewness of the tSZ fiuctuation, (r^(n)) scales 
approximately as (Tg\ whereas the amplitude of the bispectrum 
scales as cr^ with a = 11-12, as shown by Bhattacharya et al. 
(2012). In the following we do not consider the dependency of 
the bispectrum and the unnormalized skewness on other cosmo¬ 
logical parameters, since such dependencies are expected to be 
significantly lower than for erg (Bhattacharya et al., 2012). 

We derive constraints on erg by comparing the measured un¬ 
normalized skewness and bispectrum amplitudes with those ob¬ 
tained from simulations of the tSZ effect. The tSZ contribution 
was obtained from a hybrid simulation including a hydrody¬ 
namic component for z < 0.3 plus extra individual clusters at 
z > 0.3, and with erg = 0.789. This approach is strongly lim¬ 
ited by systematic uncertainties and the details of the theoretical 
modelling (see Hill & Sherwin, 2013). Uncertainties due to fore¬ 
ground contamination are computed using the simulations and 
are accounted for in the final error bars. 

We obtain erg = 0.77 for NILC and erg = 0.78 for MILCA. 
Combining the two results and considering model and fore¬ 
ground uncertainties we obtain erg = 0.78 + 0.02 (68% C.L.). 
Notice that the reported uncertainties are mainly dominated 
by foreground contamination. However the model uncertainties 
only account for the expected dependence of the unnormalized 
skewness upon erg, as shown in Appendix A. We have neglected, 
as was also the case in Wilson et al. (2012), the dependence on 
other cosmological parameters. We have also not considered any 
uncertainties coming from the combination of the hydrodynam¬ 
ical and individual cluster simulations. Because of these limita¬ 
tions, our error bars might be underestimated. 
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Fig. 19: Marginalised likelihood distribution for 
for tSZ and CMB based analyses. We represent the tSZ power 
spectrum analysis results assuming a mass bias, b, of 0.2 (red) 
and 0.4 (orange), the cluster number count analysis results 
(green; Planck Collaboration XXIV, 2015), and the combined 
Planck CMB and BAO analysis (Planck Collaboration XIII, 
2015) with (cyan) and without (blue) extra lensing constraints. 


7.3.2. Fit of the 1D PDF distribution 

We also derived constraints on erg by fitting the ID PDF ob¬ 
tained in Sect. 6.1. Here, we follow the formalism described in 
Hill et al. (2014) that evaluates the tSZ ID PDF theoretically 
integrating across individual cluster contributions. We use the 
Tinker et al. (2008) mass function and the Amaud et al. (2010) 
pressure profile. The later is normalised following the Y-M scal¬ 
ing relations described in Planck Collaboration XX (2014), and 
considering a mass bias parameter ofb = 0.2. All cosmological 
parameters are fixed to the Planck 2015 CMB analysis best-fit 
values, and we fit only for erg. 

We note that the Hill et al. (2014) formalism explicitly ne¬ 
glects any effects due to overlapping clusters along the line of 
sight. For this reason, and given the uncertainties in the mod¬ 
eling of the foreground residuals, the best-fit solution to the ob¬ 
served ID PDF is fitted in the region which is dominated by non¬ 
overlapping cluster and where the noise and foregrounds contri¬ 
butions are minimal. In our case, we use y > 4.5 x 10“^. 

The confidence limits on erg are obtained from a maximum 
likelihood approach, in which the likelihood has a multivariate 
Gaussian shape with a covariance matrix which only depends on 
(Tg. This covariance matrix is evaluated numerically, account¬ 
ing only for Poisson terms (both for pixel to pixel and due to 
the correlations introduced by the cluster’s y-profile). We obtain 
(Tg = 0.77 + 0.02 (68 % C.L.) both for NILC and MILCA maps 
including statistical and systematic uncertainties. 

7.4. Comparison to other Planck cosmological probes 

We have shown in Figure 17 that the amplitude of the tSZ power 
spectrum measured in this paper and from the Planck CMB anal¬ 
ysis are in good agreement. However, the Planck 2013 results 
(Planck Collaboration XX, 2014) have shown tension between 
CMB and tSZ derived constraints on (Tg for wide range of ex¬ 
periments including Planck . Figure 19 shows the marginalised 
likelihood distribution for (Tg(Qni/0-28)^/^ as obtained from the 
combined Planck CMB and BAO analysis (Planck Collaboration 


XIII, 2015) with (cyan) or without (blue) lensing constraints. 
We also presents the results obtained for the Planck 2015 clus¬ 
ter number count analysis assuming a mass bias of 0.2 (green; 
Planck Collaboration XXIV, 2015), and for the tSZ power spec¬ 
trum analysis in this paper assuming a mass bias of 0.2 (red) 
and 0.4 (orange). We observe that assuming a mass bias of 0.2 
the two tSZ analyses are in good agreement but the tension 
with the CMB measurements remains. This tension can be al¬ 
leviated by assuming a larger mass bias that increases the value 
of (Tg(Qm/0.28)^/^. Furthermore, we find that including lensing 
constraints leads to smaller values of (Tg(Dni/0.28)^/^, which as 
a consequence are in better agreement with the tSZ results. From 
this we can conclude that the LSS data would have a marginal 
preference for small values of the mass bias. 

8. Summary and conclusions 

Because of its wide frequency coverage from 30 to 857 GHz, 
the Planck satellite mission is particularly well suited for the 
measurement of the thermal Sunyaev-Zeldovich effect. Working 
with the Planck frequency channel maps from 30 to 857 GHz, 
we have reconstructed the tSZ signal over the full sky using 
tailored component separation methods. 

We tested and validated the Planck y-maps extensively and 
characterized them in terms of noise and foreground contamina¬ 
tion. As expected the noise in the y-map is inhomogeneous and 
can be characterized by pixel dependent variance and a homoge¬ 
neous correlated Gaussian noise. Foreground contamination by 
thermal dust emission is found to be important at large angu¬ 
lar scales. Additional foreground contamination is due to radio 
and IR point sources for which a mask is provided. In terms of 
tSZ signal we find good agreement between the fiux of blindly 
detected clusters in the y-map and that measured for clusters in 
the Planck cluster sample. Furthermore, we find that the sensi¬ 
tivity of the y-map is sufficient to detect faint and diffuse struc¬ 
tures such as bridges between merging clusters. Moreover, we 
have proved via a stacking analysis that the very low signal-to- 
noise regions in the y-map preserve the tSZ signal even for small 
galaxy groups (tens of galaxies). 

After accounting for foreground contribution, mainly ther¬ 
mal dust emission at large angular scales, and clustered CIB and 
point sources at small angular scales, we have derived from the 
y-map the tSZ angular power spectrum in the multipole range 
from 9 < i < 1411. This extends significantly the range with 
respect to previous measurements (Planck Collaboration XXI, 
2014) giving for the first time access to the 2-halo term contri¬ 
bution. The cosmological analysis of the tSZ power spectrum 
allows us to set constraints on cosmological parameters rep¬ 
resenting matter content in the Universe, mainly (Tg and 
These constraints are consistent with those obtained from clus¬ 
ter number counts (Planck Collaboration XXIV, 2015) and in 
soft tension with those derived from CMB analysis (Planck 
Collaboration XIII, 2015). 

The analysis of the non-Gaussian properties of the y-map 
using the ID PDF, the unnormalized skewness and the bispec¬ 
trum of the map have confirmed the tSZ nature of the signal. 

The Planck y-maps and additional ancillary data (noise vari¬ 
ance maps, foreground masks and ILC weights) are made avail¬ 
able to the public for the Planck 2015 release (see Appendix C 
for details). These y-maps are expected to be useful in a wide 
range of astrophysical and cosmological analyses with clusters. 
For any of these analyses, and depending on the scientific goals. 
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the inhomogeneous properties of the noise, and the systematics 
and foreground contamination should be taken into account in 
different ways as described in this paper. Regions masked by 
the point source mask should never be used. In the case of pixel 
based analyses quality flags can be defined by combining the 
information from the variance map and the various foreground 
masks. For power spectrum, cross correlation and higher order 
statistic analyses we remind the fact that the y-maps present sig¬ 
nificant foreground contamination that needs to be taken into ac¬ 
count both by masking highly contaminated regions (namely the 
Galactic plane region) and by using adequate foreground mod¬ 
els to which the ILC weights are applied. Taking these necessary 
precautions, the Planck y-maps will prove a very useful tool for 
the community. 
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The 2-halos term (Komatsu & Kitayama, 1999; Diego & 
Majumdar, 2004; Taburet et al., 201 1) is given by: 


/^2halos 
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Afmin 


dM 
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where P(k,z) is the 3D matter power spectrum at redshift z. 
B(M, z) is the time-dependent linear bias factor that relates the 
matter power spectrum, P(k,z), to the power spectrum of the 
cluster correlation function. Following Komatsu & Kitayama 
(1999, see also Mo & White 1996) we adopt B(M,z) = 1 + 
ly^(M,z) - l)/Sc(z), where v{M,z) = 6c{M)ID{z)(r{M), (t{M) is 
the present-day rms mass fiuctuation, D{z) is the linear growth 
factor, and ^c(^) is the threshold over-density of spherical col¬ 
lapse. 

Finally, we use the Tinker et al. (2008) mass function, 
dn(M, z)ldM, including an observed-to-true mass bias b, as dis¬ 
cussed in detail in Planck Collaboration XX (2014) , and we 
model the SZ Compton parameter using the pressure profile of 
Arnaud et al. (2010). 


Appendix A: Modelling the expected tSZ signal 

A.1. tSZ power spectrum 

The representation of the y-map in spherical harmonics, 
reads 

y(n) = X ytm km(«)- (A.l) 

im 

Thus, its angular power spectrum is given by 


A.2. Mh moment of the tSZ field 


Assuming a Poisson distribution (1-halo term) of the clusters 
on the sky and neglecting clustering between clusters the A^th 
moment of the tSZ signal (Komatsu & Kitayama, 1999; Wilson 
et al., 2012; Planck Collaboration XX, 2014) reads 
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where y(0, M, z) is the integrated Compton parameter along the 
line of sight for a cluster of mass M at redshift z. 


m 

Note that is a dimensionless quantity here, like y. 

As in (Planck Collaboration XX, 2014) the tSZ power spec¬ 
trum is modelled using a 2 -halo model to account both for intra¬ 
halo and inter-halo correlations: 

xCf = Cf-f Cf(A.3) 
Following (Komatsu & Seljak, 2002) the 1-halo term reads: 


A.3. Bispectrum 

The angular bispectrum is given by 

=<3'X.m.yWX3m,), (A. 8 ) 

where the angle-averaged quantity in the full-sky limit can be 
written as 

Z XkSr- <A.9) 
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where dVdidzdQ) is the comoving volume per unit redshift and 
solid angle and n(M, z)dM dVc/(dzdD) is the probability of hav¬ 
ing a galaxy cluster of mass M at a redshift z in the direction 
dkl. The quantity y^ = y^(M, z) is the 2D Fourier transform on 
the sphere of the 3D radial profile of the Compton y-parameter 
of individual clusters. 
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and satisfies the conditions mi -f m 2 + m 3 = 0 , -f ^2 + ^3 = even, 
and \^i -£j\ < £k ^ + ^ 7 , for the Wigner 3 j function in brack¬ 

ets. Assuming a Poissonian spatial distribution of the clusters as 
above, the bispectrum reads (Bhattacharya et al., 2012) 
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Appendix B: Bispectrum cosmic variance 


where v = rjr^, 4 = JJdz)lr^, is the scale radius of the 3D Following (Lacasa, 2014, chapter 2), the bispectrum 
pressure profile, Da(z) is the angular diameter distance to red- cosmic variance is composed of a Gaussian term, a 
shift z and Pq is the electron pressure profile. bispectrumxbispectrum term, a spectrumxtrispectrum term 
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and a connected 6-point term. Due to the lack of model or mea¬ 
surement of the trispectrum and 6-point function, we neglected 
the last two terms. Note that they are however expected to yield 
a subdominant contribution. Thus we have in full-sky : 


- Gaussian cosmic variance 
Q 2 Q 3 


VarG(/?^i4^3) = 




( 6 equilateral 
2 isosceles 
1 general 


where 
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An 


0 0 0 
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and where Q is the auto power spectrum of the Compton 
parameter map, thus containing the noise contribution. 

- Bispectrumxbispectrum cosmic variance 
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This is the only term which gives off-diagonal contributions 
to the covariance matrix. 


For our purpose, this cosmic variance is multiplied by /sky and 
binned to the appropriate binning scheme. 

We also consider systematic errors induced by foreground 
residuals or masking effects. We estimate systematic errors due 
to component separation uncertainties from the half difference of 
the NILC and MILCA bispectra. Masking effects are normally 
corrected for using simulations, which may under- or overesti¬ 
mate leakage from large to small scales. We thus take a conser¬ 
vative ±25% error on the debiasing ratio, consistently with the 
fact that the selected configurations have a ratio within +25% of 
/sky ^(^i)^(^ 2 )^(^ 3 )- This error is most likely a conservative 
overestimate. 


Appendix C: Products 

In the following we list the y-map related products delivered in 
the Planck 2015 data release'll 

- Full-sky MILCA and NILC y-maps for the full mission and 
for the first (F) and second (L) halves of Planck stable point¬ 
ing period in Compton parameter units (see Section 3.2). 

- ILC weights per filter and per frequency used for the recon¬ 
struction of the MILCA and NILC full mission y-maps, as 
described in Section 3.2. 

- Variance map accounting for the non homogeneous coverage 
and power spectrum of the correlated homogeneous counter¬ 
part, C^, for the MILCA and NILC full mission y-maps in 
Compton parameter units (see Sect. 4.2). 

- Point source masks including known radio and IR sources as 
described in Section 4.4.1. 

- Galactic masks used in the analyses presented and in 
Sections 5 and 6. 

^ A more detailed description is given in the Planck explanatory sup¬ 
plement. 
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